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ABSTRACT 


A  straightforward  and  general  program  of  spectral  analysis  has  been  written  in 
FORTRAN  II  for  the  IBM  7090  and  is  described  in  this  report.  The  analysis  is  valid 
for  sequences  of  equidistant  data  sampled  from  realizations  of  second  order  stationary 
stochastic  processes.  Alternatively,  the  program  may  be  used  to  estimate  the  transfer 
function  gain  characteristic  of  a  linear  system  on  the  basis  of  its  sampled  output. 

The  Inputs  to  the  program  consist  of  the  data  sequence  to  be  analyzed  and  FOUR 
control  parameters. 

The  Output  consists  of  listings  of  the  estimated  values  and  of  three  CALCOMP 
plots  of 

(i)  The  sample  auto  covariance  (ACV)  functions. 

(ii)  The  power  spectral  density  (PSD). 

(iii)  A  log-log  plot  of  the  PSD. 

The  estimated  PSD  is  consistent,  being  a  periodogram  smoothed  with  Hanning  weights. 

After  a  brief  introductory  discussion  in  Section  I,  Section  II  proceeds  with  a 
sketch  of  the  analytic  background,  and  a  discussion  of  the  parameters  critical  to  a 
power  spectral  analysis.  Section  III  is  a  description  of  the  program  and  has  some 
sample  estimates.  Operating  instructions  are  given  in  Section  IV  as  well  as  a  complete 
description  of  the  outputs.  Section  V  describes  an  alternate  use  of  the  program. 

Section  VI  points  to  possible  modifications  to  tailor  the  program  more  nearly  to  the 
individual  requirements  of  a  prospective  user. 

References  are  listed  in  the  back  of  the  report  and  are  noted  in  the  text  by 

indicating  the  reference  number  in  square  brackets.  For  someone  anxious  to  use  the 

program  it  may  be  best  to  read  the  section  on  operating  instructions  (Section  IV)  first 

along  with  Section  II  C  which  described  the  significance  of  the  program  parameters. 

Accepted  for  the  Air  Force 
Stanley  J.  Wisniewski 
Lt  Colonel,  USAF 
Chief,  Lincoln  Laboratory  Office 
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I. 


INTRODUCTION 


The  last  15  years  have  seen  a  renewed  interest  in  the  estimation  of  power  spectra 
as  a  tool  in  the  periodic  analysis  of  data.  One  of  the  primary  concerns  has  been  to  find 
suitable  modifications  of  periodogram  analysis  to  render  consistent  ( in  the  statistical 
sense)  the  estimates  of  power  spectra  (or  power  spectral  density),  and  the  principal 
device  used  to  obtain  consistency  has  been  smoothing  of  the  periodogram 3  A  multitude 
of  smoothing  sequences  (or  equivalently  lag  windows)  have  made  their  appearance  and 
are  treated  in  varying  detail  in  Refs.  [2,  3,  4] .  There  still  is  a  good  deal  of  controversy 
over  the  problem  of  choosing  the  parameters  of  the  smoothing  process  (bandwidth  or 
equivalently  the  truncation  point  in  the  corellogram)  to  insure  an  optimal  estimate,  let 
alone  over  initial  agreement  on  a  criterion  of  optimality. 

The  basic  dilemma  is  that  given  only  a  sample  of  data,  and  no  background 
information  on  either  the  process  generating  the  data  or  the  recording  device  used  in 
gathering  the  data,  there  exists  no  unambiguous  method  of  obtaining  an  estimate  of  the 
Power  Spectrum  which  may  be  regarded  as  best,  see  Ref.  [  5] . 

In  the  absence  of  clearcut  guide  lines,  we  have  settled  on  a  moderately  simple  and 
practical  procedure  ( periodogram  smoothing  with  Hanning  Weights)  and  have  incorporated 
it  into  a  program  in  such  a  way  that,  should  another  smoothing  sequence  appear  over¬ 
whelmingly  preferable,  the  modification  to  the  program  can  be  made  most  trivially. 

The  key  note  of  the  program  (the  AMPSDE  Program:  Any  Man’s  Power  Spectral 
Density  Estimator)  is  its  simplicity  and  ease  of  use.  The  number  of  available  options 
has  been  reduced  to  a  minimum  without  curtailing  the  program's  general  usefulness  as 
an  exploratory  device.  As  previously  mentioned,  and  as  described  more  fully  below, 
minor  modifications  are  easy  to  make  by  an  interested  user  in  the  light  of  his  own  personal 
requirements,  principally  in  the  areas  of  input,  output  and  choice  of  smoothing  windows. 

t  An  introductory  exposition  of  the  general  problems  encountered  in  the  analysis  of  data 
from  the  point  of  view  of  its  frequency  content  can  be  found  in  Ref.  [  1  ] . 
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In  a  similar  vein,  the  outputs  of  the  program  are  provided  without  option,  to 

insure  a  minimum  of  complexity  in  use  and  maximum  ease  in  interpretation  of  results. 

1* 

For  each  set  of  data  to  be  analyzed,  three  plots  1  are  invariably  provided: 

(i)  A  normalized  Auto-Covariance  ( ACV)  plot. 

(ii)  A  Power  Spectral  Density(PSD)  plot  on  linear  scale. 

(iii)  A  Power  Spectral  Density(PSD)  plot  on  a  log-log  scale. 

Along  with  the  graphical  outputs,  the  program  prints  a  number  of  lists  of  the  estimated 
values  which  are  more  fully  detailed  later  in  the  report. 

The  inputs  to  the  program  consist  of  four  constants,  in  addition  to  the  data 
sequence  to  be  analyzed.  Several  cases  may  be  stacked  in  a  single  run. 

The  next  section  proceeds  with  the  analytic  background  of  the  report,  starting  with 
a  brief  outline  of  the  continuous  case,  its  adaptation  to  discrete  data  and  ending  with  a 
discussion  of  the  parameters  which  were  selected  to  control  the  program  and  their 
interpretation. 

Section  III  is  a  detailed  description  of  the  program,  with  flow  charts  included 
wherever  they  clarify  the  presentation. 

Operating  instructions  are  given  in  complete  detail  in  Section  IV  and  a  complete 
sample  of  annotated  outputs  in  included. 

An  application  of  the  program  to  the  estimation  of  transfer  functions  is  discussed 
in  Section  V. 

Finally  Section  VI  points  out  the  ease  with  which  the  program  may  be  modified. 
Appendix  A  contains  a  brief  note  on  discrete  Fourier  Cosine  Transforms.  A 
complete  program  listing  is  given  in  Appendix  B.  Appendix  C  exhibits  a  sample  of 
listed  outputs,  while  graphical  outputs  are  illustrated  in  Figs.  3,4  and  5. 

II.  ANALYTIC  BACKGROUND 

The  estimation  procedure  used  in  the  program  follows  the  simplest  lines 
suggested  in  Ref.  [3] ,  but  there,  this  simplicity  is  somehow  imbedded  in  such  a  wealth 
of  information  that  the  clarity  of  the  procedure  is  all  but  apparent. 

t  The  graphical  outputs  are  produced  on  a  CALCOMP  plotter. 
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The  basic  theorem  states  that  if  a  process  x^  has  a  summable  ^  auto  covariance 
sequence  y(k),  k  =  0,  ±1,  ±2,  . . . ,  then  its  power  spectral  density  exists  and  is  given  by 

OO 

=  2?  X  Y(k)  cos(ka,)* 

k=-°° 

All  estimates  f(co)  of  f(w)  which  have  been  devised  so  far  are  ultimately  of  the  form 

N-l 

f(w)=  27  Yj  hT(k)<^T(k)cos(koj), 
k=-N+l 

where  ^^(k)  is  an  estimate  of  y(k)  based  on  the  sample  (dependence  on  sample  size  has 
been  shown  by  the  subscripted  T  =  NAt)  where  n  =  0, 1,  2, . . . ,  N  and  is  generally  taken 
as: 


<PT(k)  = 


1 

N  -  |"k  | 


N—  |  k  | 

Z  <Xn~*><xn+|  k  | 


n=0 


k  =  0,  ±1.  ±2 . N-l 


V  0  ,  |k|>N 

and  where  h^(k)  is  some  convariance  averaging  kernel  depending  on  the  lag  index  k 
and  possibly  a  sample  size  as  well.  Many  different  forms  of  {h^(k)}  have  been  suggested, 
but  no  dearcut  method  exists  for  choosing  an  optimal  sequence  in  the  absence  of  addi¬ 
tional  information  on  either  the  process  or  the  sampling  procedure. 


t  Summable  y(k)  means  £  I  yOO  |  < 00  • 

k=0 
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Nor  can  an  optimality  criterion  be  chosen  in  vacuo,  see  Ref.  [6];  for  example,  if 
we  wanted  to  estimate  the  PSD  at  a  single  point  we  might  seek  to  minimize 

(f[£(w0)  -f(“0)]2  . 

If,  on  the  other  hand,  one  is  interested  in  estimating  the  spectrum  over  an  interval 
a  <co  <b,  one  might  seek  to  minimize 

J  (f  |f  -  f  |  2 doj  . 
a 

It  is  perhaps  simplest  to  begin  with  a  brief  exposition  of  the  continuous  case  and  proceed 
to  specialize  to  the  situation  which  generally  arises  in  practice: 

A.  Continuous  Case 

The  simplest  situation  obtains  when  x(t)  is  a  zero  mean,  stationary,  ergodic, 
stochastic  process  whose  auto  covariance  function  (ACV) 

T/2 

(1)  t(t)  =  lim  -  C  x(t)x(t+r)  dt  =  yf-r) 

T— -°°  J 

-T/2 

admits  a  spectral  representation  in  terms  of  f(«),  where 

oo 

(2)  f(w)  e  1WT  t )dr 


f(w)  is  called  the  power  spectral  density  of  the  process  and  can  be  written  as 

T/2 


(3)  f(w)  =  lim  - 

T— *°o  1 


I 


x(t)e  ldt 


-T/2 


Equations  (1)  and  (2)  express  the  Wiener-Khinchin  theorem,  while  Eq.  (3)  shows  why  the 
periodogram  apprears  to  be  a  natural  estimate  of  the  PSD,  (Ref.  [  1]  ). 
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It  may  be  worth  placing  the  ACV  function  and  the  PSD  of  x(t)  in  evidence  as  a 
Fourier  pair: 


Y(t)  =  2^  y  f(w)e1UjTdcu 


(4) 


y(T)e  dT 


I 

-oo 

Since  both  y(-r)  and  f(w)  are  even  functions,  we  may  rewrite  (4)  as 


(5) 


>,<T)  =  -h  I  f(w)  COS(WT)dW  =  _1_  y  f(w)  COS(WT)dW 


f(w)  =  ^  y(-r)cos(a>)dT  =  2  ^  y(r)  cos(a>-r)dT  , 
-°o  0 


where  the  second  equality  expresses  y(T)  and  f(o>)  as  one-sided  integrals. 

B.  Discrete  Case 

In  order  to  apply  digital  methods  we  must  sample  x(t)  at  equidistant  intervals 
At  over  a  finite  length  of  time  T,  to  obtain  a  sequence 


where 


Vxi . Xn”,,,XN 

x  =  x(nAt)  and  N  =  T/At. 


As  a  first  step,  in  practice,  it  is  necessary  to  remove  the  DC  component  (nonzero 
mean)  from  the  data  before  computing  the  sample  autocovariance,  C^.  This  is  readily 
accomplished  by  replacing  the  original  {x^}  by  {x^  -  x},  where 
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N 


n=0 


One  proceeds  next  to  estimate  the  process  auto  covariance  by  the  sample  auto  covariances. 


(6) 


V 


N-|k| 


N  -  ^0 


n+  |M 


|k|  =  0,1, ...,K 


where  K  is  the  maximum  lag  index  (necessarily  <N). 

The  obtained  sequence  of  sample  ACV's  is  normalized  to 

C(0)  =  C(1 )  C(K ) 

0  C(0)  ’  1"C(0)  ’•••’  K“C(0)  , 

and  a  discrete  finite  cosine  transform  of  the  resulting  sequence  is  taken,  as  shown  in 
the  following  relations 


K-l 

p<rAt<co  +  2I  Ck+CK} 

k=l 

K-l 

P,  =  At{C_  +  2  /  C.  cos  —  +  CT,  cos  tt} 

1  1  0  Li  k  K  K 

k=  1 

+  cos^+CK-2’) 

k=  1 

•  •  •  • 

(7)  ... 

•  •  •  • 

K-l 

P  =  AtjC,.  +  2  /  C,  cos  +  CT,  cos  rtf  } 

r  0  Li  k  K  K 

k=l 


6 


K-l 


P.,  =  At{C_  +  2  /  C  cos  kTr  +  C  cos  Kir} 
K  U  Lj  K  K 


k=l 


Finally,  the  results  may  be  smoothed  by  Hanning  weights 


U0  =  2<P0+P1> 


<8)  Uk  =  ?  Pk-1  +T Pk  +  i  Pk+1 


UK=2<PK-1+PK> 


k  =  1,  2, . . . ,  K-l 


C.  Choice  of  Parameters 

A  cursory  examination  of  the  previous  section  (II,  B)  reveals  that  three 
parameters  govern  the  process  of  generating  PSD  estimates.  These  are: 

(i)  The  number  of  data  points  in  the  sample  (N+l). 

(ii)  The  sampling  rate  or  the  inter-sample  time  interval,  (DT). 

(iii)  The  maximum  lag  index  for  calculation  of  auto  covariances,  (K). 

For  purposes  of  plotting  and  improved  resolution  in  the  vertical  scale  a  fourth  parameter 
has  been  included  in  the  program: 

(iv)  A  scale  factor  by  which  all  PSD  estimates  are  multiplied. 

We  will  discuss  the  significance  of  the  scale  factor  in  Section  III  where  the  calculations 
and  output  scaling  performed  by  the  program  are  described.  We  now  turn  to  a  discussion 
of  the  parameters  which  have  analytical  significance. 

The  most  important  single  parameter  of  the  program  is  the  sampling  rate  as  it 
controls  the  highest  frequency  component  of  the  data  which  can  be  meaningfully  identified. 
The  maximum  lag  index  determines  the  resolution  in  PSD  obtainable  by  the  program. 
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The  number  of  data  points  determines  the  quality  of  the  estimates,  depending  on  the 
signal -to -noise  ratio  in  the  data  and  the  data  collection  scheme. 

1.  DT  and  the  Sampling  Theorem  -  The  sampling  theorem  states,  in  one  of  its 
forms,  that  if  At  seconds  is  the  sampling  interval  in  a  sequence  of  data,  then  the 
frequency  components  of  the  process  generating  the  data  with  frequency  w  >  2. 

7 r 

radians  /second  cannot  be  distinguished  from  these  with  frequency  in  the  range  (0,  — ) 

on  the  sole  basis  of  the  sampled  values,  see  Ref.[7], 

As  a  consequence  it  is  best  to  choose  At  small  enough,  in  practical  situations,  to 

insure  that  negligible  power  is  contained  by  the  process  beyond  the  frequency 
7 r 

WNy  =  At  ^Ny  =  fre<luency)*  A  high  sampling  rate,  however,  (small  At )  means 

large  quantities  of  data;  to  avoid  unmanageable  quantities  of  data  it  may  be  worthwhile 
passing  the  process  through  a  low  pass  filter  before  sampling,  if  it  is  felt  that  there 
exist  high  frequency  components  in  the  process  which  are  of  no  interest  to  the  analyst. 

In  most  situations,  however,  the  analyst  is  presented  with  a  collection  of  data 
where  the  sampling  interval  has  already  been  determined.  In  this  case,  one  must  set 
the  parameter  DT  equal  to  At  and  there  is  no  way  of  determining  whether  there  has  been 
any  aliasing  of  high  frequencies,  see  Ref.  [  1] . 

7 T 

The  program  automatically  restricts  itself  to  the  basic  interval  from  0  to  — 
in  the  calculation  of  estimated  PSD. 

2.  K  and  PSD  Resolution  -  Equation  (7)  of  Section  II.  B  can  be  rewritten  in  matrix 
form  as 

P  =  M  C 

where  P  is  a  (K+l)  vector  of  estimate  PSD  values,  M  is  the  coefficient  matrix  of  cosines, 
and  C  is  the  (K+l)  vector  of  normalized  sample  ACV’s. 

These  K+l  values  of  PSD  are  equally  spaced  over  the  interval  [  0,co^^) ,  so  that  if 
Aw  is  the  angular  frequency  resolution  desired,  the  analyst  must  set  K  such  that 
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_N> 

K 


Aw 


f 


Ny  ~  Af 
K 
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or 

/q\  K  ~  =  i  .  K  ~  V  = _ L_ 

K  '  Aw  Aw  (At)  ’  Af  2Af(At) 

3.  Sample  size  and  S/N  Ratio  -  The  accuracy  of  estimation  is  directly  dependent 
on  sample  size.  If  the  signal -to -noise  ratio  does  not  render  the  sampled  values  useless, 
a  larger  sample  will  in  general  yield  better  estimates.  It  is  also  advisable  not  to  run 
the  maximum  lag  at  which  ACV's  are  calculated  too  close  to  the  sample  size;  indeed, 
Eq.(6)  shows  that  as  the  lag  index  grows,  the  number  of  terms  in  the  estimate  of  C, 
becomes  smaller  and  hence  the  sample  ACV’s  become  poorer  estimates  of  the  process 
autocovariance  for  large  lags.  As  a  rule  of  thumb,  a  reasonable  choice  for  maximum 
lag  index  might  be  20%  of  the  sample,  so  that,  if  K  is  determined  on  the  basis  of 
resolution  requirements,  one  would  let 

N  *  5K. 

4.  Summary  -  In  recapitulation,  after  the  sampling  interval  DT  is  set,  the 
maximum  frequency  content  of  the  sample  is  determined  by 

7T  1 

w  =  —  radians/ second  or  f  =  cycles/second 

The  analyst  then,  as  a  first  approximation,  may  choose  a  resolution  level  by  setting 
the  maximum  lag  index  K  according  to  Eq.  (9),  and  the  sample  size  (N+l)  to  exceed  5 
times  the  maximum  lag.  For  high  S/N  ratios,  K  may  go  as  high  as  50%  of  sample 
length.  In  threshold  situations,  one  may  set  K  to  be  only  5%  of  the  sample  length  in 
order  to  stabilize  the  autocovariance  estimate. 

III.  PROGRAM  DESCRIPTION 

A  flow  chart  of  the  program  is  shown  in  Fig.  1.  A  number  in  parentheses  next  to 
a  box  refers  to  the  equation  in  the  text  relevant  to  the  process  described  in  the  box.  An 
encircled  number  corresponds  to  the  output  described  in  Section  IV.  D.  Flow  of  control 
is  marked  in  solid  lines;  flow  of  data  is  shown  in  dotted  lines. 
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Fig.  1  Program  Flow  Chart 


The  program  is  written  in  FORTRAN  II  and  makes  use,  for  its  graphical  outputs, 
of  a  set  of  CALCOMP  plotting  subroutines  which  are  documented  in  Ref.  [8] .  Apart 
from  these  subroutines,  the  program  is  completely  straightforward  as  can  be  seen  from 
the  flow  chart,  or  by  consulting  the  listing  given  in  Appendix  B. 

As  has  already  been  noted,  the  program  is  designed  to  spare  the  user  from  the 
necessity  of  choosing  one  out  of  an  elaborate  list  of  options.  Instead,  he  will  have  to 
make  judicious  choices  of  the  four  parameters  controlling  the  program. 

We  will  now  proceed  to  describe  briefly  several  aspects  of  the  program  which 
warrant  amplification. 

A.  Inputs 

Four  constants  completely  determined  the  course  of  the  program.  These 

are 

(i)  The  maximum  sample  index  N  equal  to  one  less  than  sample  size, 

(ii)  The  maximum  lag  index  K  for  which  sample  autocovariances  are 
calculated, 

(iii)  The  sampling  interval  DT, 

(iv)  A  scaling  factor  SF. 

The  significance  of  these  parameters  (except  for  the  4th)  has  been  examined  in 
Section  II.  C. 

B.  Computations 

The  computations  follow  closely  the  procedure  outlined  in  Section  II.  B.  The 
dc.  component  is  removed  from  the  data  before  the  frequency  content  of  the  data  is 
examined.  Autocovariances  are  calculated  in  the  standard  fashion.  The  cosine 
transform  of  the  ACV  function  is  taken  and  then  smoothed,  and  final  outputs  are  made 
while  auxiliary  outputs  are  generated  at  various  points  of  the  program. 

A  word  is  in  order  concerning  the  program's  handling  of  Eq.(7),  which  can  be 
rewritten  as 

P  =  MC 
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where  M  is  a  (K+l)x(K+l)  matrix  of  cosines.  A  ingenious  scheme  avoids  the  necessity 

9 

of  calculating  and  storing  (K+l)  elements.  Instead,  only  K+l  cosine  values  need  be 
calculated.  Indeed,  all  the  entries  of  matrix  M  are  found  among  the  K+l  values 

k7T 

m  =  cos  ( -rr  )  ,  k  =  0,  1, 2, . . . ,  K. 

K  Jx 

These  are  the  only  values  computed  by  the  program.  The  rows  of  matrix  M  are 
generated  one  by  one  by  selecting  the  appropriate  entries  in  the  "cosine  table, "  as 
described  below.  The  cosine  table  is  a  sequency  of  2K+1  elements  formed  by  continuing 
the  sequence  symmetrically,  by  reflection  as  shown  in  Fig.  2,  for  K  =  8,  by  way  of 
illustration. 


The  rows  of  matrix  M  are  formed  by  continuing  periodically  the  sequence  stored 
in  the  cosine  table,  or  equivalently,  by  assuming  that  the  cosine  table  is  circularly 
stored  so  that  the  value  following  the  last  entry  is  the  value  stored  in  the  first  entry  - 
and  by  selecting  elements  of  M  by  the  following  scheme: 

The  first  row  of  M  contains  K+l  ones. 

The  second  row  of  M  contains  the  first  K+l  elements  of  the  cosine  table. 
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The  third  row  of  M  contains  the  every  second  element  of  the  cosine  table, 
until  a  row  or  length  K+l  has  been  completed. 

The  fourth  row  of  M  contains  every  third  element  (Mod  2K+1,  the  length 
of  the  table)  of  the  cosine  table,  until  a  row  of  length  K+l  has  been  generated. 
And  so  on,  until  (K+l)  cosine  vectors  have  been  generated  and  for  each,  the 
inner  product  has  been  formed  with  the  autocovariance  vector  (sequence) 

C,  the  resulting  (K+l)  vector  of  inner  products  being  the  periodogram. 

The  PSD  estimated  by  the  program  being  a  smoothed  periodogram,  a  smoothing  process 
occurs  next,  governed  by  Eq.  (8). 

Since  the  smoothed  estimate  has  been  calculated  from  the  normalized  auto¬ 
covariances,  the  of  Eq.(8)  are  estimates  of  a  normalized  spectral  density  P^,(f). 

This  may  be  what  is  desired  in  that  comparisons  are  sometimes  best  made  upon  a 
normalized  spectrum,  say 

(10a>  p‘k  -  W  *  KUk 

where  the  factor  k  is  independent  of  k  or  of  the  sampled  process,  but  may  depend  upon 

the  possible  overall  frequency  content  (  i.e. ,  on  At).  Many  times  one  will  desire  true 
levels  given  by 

(iob)  frp«k)-Vk>  =  '\ 

2 

where  cr  is  the  sample  variance.  Finally,  one  may  desire  true  levels  for  the  process 
but  must  introduce  some  scaling  constant  due  to  the  instrumentation  and/or  pickup 
sensors.  For  this  case,  one  desires 

(10c)  ^k=a<r2uk 

2 

where  as  before  cr  is  the  sample  variance  and  a  is  some  positive  constant  determined 
independently  of  the  program. 
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The  AMPSDE  program  is  general  enough  to  include  all  three  of  these  options  with 
the  specification  of  a  single  parameter  called  the  scaling  factor  (SF).  Thus, 

(i)  if  SF  =  0.  0,  P  =  kU  where  k  -  1/At 

K  K 

(ii)  if  SF  =  1.  0>P,  =  <r  2U  where  c r  2  =  sample  variance 

K  K 

o  o 

(iii)  if  SF  =  a>  0,  P^  =  where  again  cr  =  sample  variance 

The  particular  choice  of  k  in  the  first  (normalized)  case  has  been  made  so  that  the 
log-log  plot  of  the  spectra  will  always  have  sufficient  detail  for  all  ranges  of  frequency 
content.  Also,  if  one  desires  the  physical  (one  sided)  power  spectrum,  one  must  use 
a  scale  factor  of  two  (2.0)  in  the  second  case  or  2a;  in  the  third  case. 

The  estimated  power  levels  usually  have  a  wide  range  of  variation;  therefore  the 
program  provides  for  the  calculation  and  plotting  of  the  spectral  estimate  on  a  logarith¬ 
mic  (db. )  scale.  These  db.  levels  being  calculated  from  the  equation 

(11)  P.(fk)  =  iO  logi0(Pk)  (indb.) 

where  the  P,  are  the  scaled  smoothed  estimates.  The  smoothed  values  which  are 
k 

negative  or  which  yield  power  levels  below  -35  db.  are  set  by  the  program  to  yield 
exactly  -35db.  The  resulting  db.  power  levels  are  also  plotted  on  a  logarithmic 
frequency  scale,  the  plot  having  an  overall  range  of  60  db.  (from  -35  to  +25).  The 
scale  factor  (SF)  described  above  may  also  be  used  to  obtain  a  better  position  of  the 
spectral  signature  on  the  db.  plot.  The  utility  of  this  log-log  plot  in  system  studies  is 
discussed  in  Section  V. 

The  remainder  of  the  program  is  taken  up  by  plot  generation  and  output  processing. 
C.  Sample  Outputs 

On  the  following  three  pages,  we  give  sample  output  plots  for  three  important 
classes  of  inputs.  The  actual  input  data  had  been  previously  synthesized  from  computer 
programs . 
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N  ■  2099  DT- 0.005 

K  »  200  SF  »  2.0 


cn 


Fig.  3  Estimates  from  Sine  Wave 


C 28 -129 


N  »  2099  DT  »  0. 005 

K  «  200  SF  «  200.  0 


Fig.  4  Estimates  from  White  Noise 


C28-132 


N  =  2099 
K  *  200 


DT  *  0.  005 
SF  =0.0 


10-00  190.00 


Fig.  5  Estimates  of  Colored  Noise 


For  Fig.  3,  the  input  was  a  23  cycle,  unit  amplitude  sine  wave. 

For  Fig.  4,  the  input  was  a  set  of  Gaussian  distributed  random  numbers  of 
approximately  unit  variance.  These  represent  a  sampled  version  of  a  white  Gaussian 
noise  procee.  With  SF  =  1/DT,  the  mean  value  of  the  spectral  density  is  the  zero  db 
line. 

For  Fig.  5,  the  input  sequence  was  the  above  noise  sample  but  first  shaped  by  a 
low  pass  filter.  This  sample  is  quite  characteristic  of  real  physical  (colored)  noise. 

IV.  OPERATING  INSTRUCTIONS 

The  AMPSDE  program  is  written  in  FORTRAN  II  for  the  IBM  7090  computer.  To 
use  the  program,  one  needs: 

a)  A  program  deck  with  CALCOMP  subroutines 

b)  One  or  more  input  data  decks 

c)  FORTRAN  system  tape 

d)  Output  tapes 

et )  CALCOMP  plotter  (off-line). 

A.  Usage 

A  usage  diagram  is  shown  in  Fig.  6.  As  illustrated  there,  it  is  desirable 
to  prestore  a  run  prior  to  operation  as  it  may  be  time  consuming  to  read  in  large 
quantities  of  data  on  line.  The  output  of  a  run  consists  of  a  listing,  on  the  system 
output  tape  (A3),  and  of  a  BCD  tape  of  graphical  data  for  the  CALCOMP  plotter,  on 
tape  (A6).  Tape  (A6)  contains  one  file  per  data  case;  each  file  on  A6  contains  three 
plots.  When  requesting  CALCOMP  outputs,  it  is  necessary  to  remember  how  many 
cases  were  run,  as  each  case  occupies  one  file. 

B.  Deck 

The  composition  of  an  AMPSDE  deck,  ready  for  operation,  is  illustrated 
in  Fig.  7.  After  the  usual  FORTRAN  System  Cards,  a  binary  deck  consisting  of 
AMPSDE  and  the  CALCOMP  Subroutines  follows:  A  FORTRAN  System  — *  DATA — 
card 

t  The  program  can  be  used  without  a  plotter,  but  a  buffer  tape  A6  must  then  be  used. 
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Fig.  6  Usage  Diagram 
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CONTROL  CARD  (4  parameters ) 


*DATA 


^CALCOMP  SUBROUTINES 


'AMPSDE  PROGRAM  BINARY  DECK 


*  XEQ 


*1 D  CARD 


C28-220 


Fig.  7  Deck  Composition 


precedes  the  data  which  may  consist  of  several  cases.  Each  case  is  made  up  of  one 
control  card  and  a  deck  of  data  cards.  The  last  card  of  the  composite  deck  should 
contain  a  negative  fixed  point  integer  in  the  first  field  (110),  signaling  to  the  main 
program  that  all  cases  have  been  completed. 

C.  Card  Formats 

There  are  three  types  of  cards  to  describe.  These  are: 

(i)  Control  Cards:  This  card  contains  the  four  control  parameters  of 
the  case  in  the  format: 

Format  110  110  F  10.8  F10.5 

Variable  N  K  DT  SF 

N+l  is  the  number  of  data  points 

K  is  the  maximum  lag  index  in  calculating  autocovariances 
DT  is  the  data  spacing 

SF  is  a  scale  factor  governing  the  scale  of  the  plots. 

(ii)  Data  Card:  Each  card  contains  7  floating  point  values  in  Column  1-70 
format  F  10.  7. 

(iii)  End  of  Run  Card:  This  card  signals  that  all  cases  in  a  given  run 
have  been  processed.  It  contains  a  negative  fixed  point  integer  in  the 
first  I  10  field. 

D.  Outputs 

Sample  output  plots  have  been  described  in  Section  III.  C.  For  completeness 
we  shall  give  a  full  list  of  the  programs  output.  The  numbering  corresponds  to  the 
circled  numbers  on  the  flow  chart  in  Fig.  1.  The  first  six  outputs  are  printed  outputs. 

The  last  three  outputs  are  graphical. 

(i)  The  four  control  parameters  are  printed  as 

N  = _ 

K  = _ 

DT  = _ 

SF  = 
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(ii)  The  first  200  values  of  input  data  are  listed  sequentially,  8  per  line. 

(iii)  The  mean  and  variance  of  the  data  are  given  in  the  following  form: 

THE  MEAN  OR  DC  COMPONENT  OF  THE  INPUT  =  _ 

THE  VARIANCE  OF  THE  INPUT  = _ 

(iv)  The  normalized  sample  covariances  are  listed  8  per  line,  starting 
with  the  normalized  variance  (autocorrelation  with  lag  0),  which  is 
identically  equal  to  1. 

(v)  The  estimated  spectral  density  is  listed,  8  values  per  line. 

(vi)  The  power  levels  in  db.  corresponding  to  the  estimated  spectral 
density  are  listed.  These  are  the  values  plotted  on  the  log-log 
plot. 

(vii)  A  plot  of  the  normalized  sample  autocovariances  (autocorrelations). 

(viii)  A  plot  of  the  estimated  Power  Spectral  Density  (PSD). 

(ix)  A  log-log  plot  of  the  PSD. 

A  complete  sample  printed  output  (i  to  vi)  is  included  as  Appendix  C. 

E.  Limitations 

The  maximum  number  of  data  points  that  the  program  can  handle  is 
10,  000.  The  maximum  lag  index  that  the  program  can  accept  is  1000.  There  are  no 
other  restrictions  apart  from  those  imposed  by  Format  Statements,  e.g. ,  input  data  is 
assumed  to  be  of  the  form  F  10.  7.  A  trivial  modification  will  permit  values  larger 
than  1000  for  K  to  be  accepted  by  the  program. 

V.  APPLICATION:  ESTIMATION  OF  TRANSFER  FUNCTION  GAIN  CHARACTERISTICS 
In  this  section,  we  describe  the  use  of  the  AMPSDE  program  in  the  estimation  of 
the  transfer  function  of  both  real  and  simulated  systems.  When  we  say  transfer  function, 
a  linear  system  is  usually  implied,  but  the  program  may  be  used  just  as  readily  to 
measure  the  spectral  transfer  characteristics  of  nonlinear  systems. 
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A.  Linear  Systems 

As  is  well-known,  a  single-input/single-output  linear  system  is  completely 
characterized  by  its  impulse  response  or  its  transfer  function  H(jw).  In  general  H(jco) 
is  a  complex  valued  function  and  thus  it  is  often  convenient  to  present  the  gain  (i.e. 
modulus  |  H(jco) | )  and  phase  ((f)  =  arg  H(jco)  separately  as  functions  of  frequency.  It  is 


y(t) 


customary  in  engineering  practice  to  present  these  plots  in  a  nonparametric  form  called 
the  Bode  diagram  (Ref.[9])  wherein  the  gain  characteristic  is  plotted  in  decibels  (db. ) 
versus  a  logarithmic  frequency  scale,  and  the  phase  is  plotted  in  degrees  versus  log 
frequency. 

In  many  applications  (e.  g. ,  radar  and  sonar  systems)  the  phase  (f>  is  not  important 
and  thus  the  principal  concern  is  the  gain  characteristic  |H(j$)| .  It  is  to  the  determination 
of  this  characteristic  which  the  present  program  is  particularly  well  suited.  If,  however, 
the  phase  information  is  also  required,  it  may  be  estimated  independently  via  cross - 
correlation  techniques. 

For  the  application  at  hand,  the  utility  of  the  AMPSDE  program  hinges  upon  the 
relation  between  the  inout/ output  spectral  densities  for  linear  systems;  namely 

S  (f)=  | H(jw)| 2  S  (f)  (0  =  2rci 

y  x 

where  S  (f)  is  the  spectral  density  of  the  input  and  S  (f)  is  the  output  spectral  density, 
x  y 

Solving  for  log  |  H(jw  | ,  we  have 

2  1og|H(j«>)|  =  logSy(f)  -  logSx(f) 


t  Cf.  Ref.  [10],  Eq.  9-36 
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or  better,  in  db's. 


20  log  |H(j«)|  =  10  log  S  (f)  -  10  log  Sx(f)  , 

and  thus  we  may  obtain  the  gain  characteristic  as  the  difference  of  two  log  spectral 
estimates. 

In  the  application  of  the  program  to  the  determination  of  Bode  gain  character¬ 
istic  for  a  real  physical  system,  one  excites  the  system  with  some  broad-band  excitation 
(noise  or  possible  an  impulse)  and  records  both  the  input  and  output  of  the  system.  The 
data  may  be  taken  directly  in  digital  form  or  as  analogue  tape  recordings  which  are 
subsequently  converted  to  digital  form  by  an  analog-to-digital  (A/D)  converter  system. 

In  either  case,  the  program  is  then  used  to  estimate  the  spectra  of  both  the  input  and 
output  record  and  the  gain  charcteristic  is  obtained  by  differencing  the  log-log  estimates. 
This  final  differencing  is  not  done  by  the  program  as  it  might  not  be  necessary  to 
calculate  both  estimates.  This  is  especially  true  for  those  applications  in  which  the 
excitation  signal  is  preserved  and  reused  in  many  spectral  determinations.  Moreover, 
(cf. ,  next  paragraph),  since  we  are  free  to  choose  the  input  x(t)  in  the  study  of  simulated 
systems,  it  may  well  be  that  Sx(f)  is  known  and  need  not  be  estimated. 

One  of  the  principal  applications  of  the  program  to  date  by  the  authors  has  been 
the  analysis  of  simulated  systems  and  the  evaluation  of  various  digital  filtering  schemes. 
In  these  applications,  the  simulated  system  or  filter  is  excited  by  a  known  sequence 
(deterministic  or  random)  for  which  the  spectrum  S^(f)  is  known  or  precalculated.  Two 
very  important  excitation  functions  in  system  studies  are  the  inpulse  function  and  white 
noise,  both  of  which  have  flat  spectra  (i.e. ,  S^(f)  =  constant).  In  these  cases,  the  Bode 
plot  is  determined  to  within  a  factor  by  the  single  estimation  of  the  simulation's  output 
spectra.  The  program  has  been  very  successfully  used  in  this  capacity  to  evaluate  some 
simple,  yet,  very  efficient  recursive  schemes  for  digital  filtering.  The  real  advantage 
of  the  technique  is  that  it  allows  one  to  measure  the  numerical  (e.g.,  round-off)  errors 
and  the  stability  of  a  filtering  scheme  for  which  a  direct  error  and  stability  analysis  may 
be  very  difficult. 
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As  an  example  of  the  application  of  the  program  to  the  estimation  of  Bode  gain 
characteristic,  consider  the  simple  three  element  linear  recursive  scheme 

y(n)  =  Ax(n)  +  By(n-l)  +  Cy(n-2) 

which  is  a  discrete  version  corresponding  to  the  continuous  2nd  order  system 

H<S>  =  - - ^ - T  • 

s  +  2£wqs  + 

For  the  particular  choice  of 

A  =  0. 00872  B  =  1. 97330  C  =  -  0. 98202 

which  is  the  sampled  continuous  system  corresponds  to 

cUq  =  67 r  £  =  0. 05  At  =  0. 005 

the  gain  characteristic  has  been  calculated  from  both  the  impulse  and  noise  response. 

The  results  are  presented  in  Figs.  8  and  9.  Figure  10  gives  the  spectrum  of  the  input 
noise  sample  used  to  drive  the  recursive  scheme.  The  noise  sample  was  the  same  one 
as  was  used  in  the  estimate  given  in  Fig.  4,  but  here  the  resolution  in  frequency  is 
doubled  (K  =  400  as  opposed  to  K  =  200).  By  rights,  the  estimated  gain  characteristic 
is  the  difference  of  the  curves  of  Figs.  9  and  10.  However,  since  the  overall  input 
spectrum  is  flat,  the  small  local  variations  in  the  output  spectrum  (Fig.  9)  do  not 
detract  much  from  a  meaningful  understanding  of  the  gain  characteristic. 

B.  Nonlinear  Systems 

For  nonlinear  systems,  the  principle  of  superposition  does  not  hold  and 
the  spectral  decompositions  of  a  system's  input  and  output  are  not  in  one-to-one 
correspondence.  However,  many  system  studies  do  require  measures  of  signal  distortion 
and  other  spectral  changes  made  by  a  system.  The  present  AMPSDE  program  is  also 
an  efficient  tool  in  this  area.  As  an  example,  we  give  in  Fig.  11  the  output  spectrum  of 
a  half-wave  rectifier  under  a  ten  cycle  sine  wave  excitation.  We  leave  it  to  the  reader 
to  verify  for  himself  that  the  resulting  spectral  components  are  the  correct  one  and  have 
the  correct  relative  magnitude. 
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Fig.  8  Frequency  Response  of  Difference  Equation  from  a  Unit  Impulse 
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Fig.  9  Frequency  Response  of  Difference  Equation  from  White  Gaussian  Noise 
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Fig.  10  Spectra  of  White  Gaussian  Noise 
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Fig.  11  Output  Spectra  of  Half-wave  Rectifier  Under  Sine  Wave  Excitation 


VI.  GUIDELINES  TO  POSSIBLE  MODIFICATIONS 


The  four  principal  areas  in  which  modifications  of  the  program  might  be  indi¬ 
cated  are: 

-  Input 

-  Autocovariance  estimation 

-  Smoothing  of  the  periodogram 

-  Output 


A.  Input 

It  has  already  been  noted  that  input  is  made  via  punched  cards,  and  that  prestoring 
a  run  obviates  the  objection  of  slow  read  in.  If,  however,  data  is  available  on  tape  in 
either  BCD  or  binary  form,  it  is  easy  to  modify  the  input  portion  of  the  program  to  handle 
the  new  medium. 

B.  Calculation  of  Sample  Autocovariance 

The  program  estimates  the  autocovariances  of  the  process  by  the  sample  covari¬ 
ances  using  the  formula 


C  = 


JL _ 

N-|k| 


*-\ 

I 


XX  I,  | 

n  n+|k| 


k  =  0,  1, . . . ,  K. 


Several  authors  prefer  to  use 


N- 1  k  | 

Ck  =  N  Z  XnXn+|k| 
n=0 


Other  modifications  are  possible,  amounting  to  weighting  suitably  the  autocovariances, 
before  taking  the  cosine  transform  (cf.  Ref.  [  11] ,  p.  59  ff.).  If  one  or  another  of  the 
covariance  estimates  appears  preferable  to  the  estimate  used  in  the  program,  it  is  easy 
to  replace  the  section  of  the  program  headed  "CALCULATE  AUTOCORRELATION 
FUNCTION"  with  an  alternative. 
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C.  Periodogram  Smoothing 

The  program  section  which  performs  the  smoothing  of  the  periodogram 
consists  of  precisely  five  statements.  These  can  readily  be  replaced  by  another 
smoothing  sequence.  Note  that  to  each  periodogram  smoothing  sequence,  there 
corresponds  an  autocorrelation  weighting  scheme  and  vice  versa.  It  appears,  however, 
simpler  in  digital  processing  to  do  the  smoothing  after  the  cosine  transform  has  been 
taken,  i.e. ,  in  the  frequency  domain,  rather  than  in  the  lag  domain. 

D.  Output 

At  user's  option,  minor  modifications  will  insure  the  omission  of  any  of 
the  presently  provided  outputs.  Additional  outputs,  however,  will  entail  more  work,  as 
required. 
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FINITE  FOURIER  COSINE  TRANSFORMS 
The  Fourier  Cosine  Transform  of  a  function  f(t)  is  defined  as 

oo 

(Al)  A(w)  =  I  f(t)  cos  u>t  dt. 

0 

The  corresponding  inversion  formula 

OO 

(A2)  f(t)  =  i  y  A(a>)  cos  cot  dw 

0 

also  holds.  If,  instead  of  being  defined  on  the  whole  axis,  f(t)  is  given  at  only  a  finite 
number  of  equi -spaced  points 

(A3)  fQ  =  f(0),  il  =  f(At),  f2  =  f(2At), . . . ,  fN  =  f(NAt), 

the  Fourier  transform  can  be  naturally  approximated  by  a  numerical  integration  scheme 
such  as  the  trapezoidal  rule.  Equation  (Al)  then  becomes 

N-l 

(A4)  A(w)={lfo  +  £  f  cos(jwAt)  +  j  f^c°s  NcoAt  At. 

3=1 

7 r 

According  to  the  sampling  theorem,  the  meaningful  range  of  to  is  the  interval  (0,  -^— ), 
and  since  N  +  l  values  of  f  are  available  we  can  select  N  +  l  values  of  to  according  to 

(A5)  »  k  —  0, 1,  2, ... ,  N. 

The  discrete  analog  of  (A4)  now  becomes 

N-l 

(A6)  \  =  A(  =  f0  +  fj  C°S  (  ^N~  *  +  2  *n}  At 

3=1 
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Note  that  an  ordinary  Fourier  cosine  expansion  of  the  sample  {T}  yields 


N-l 


2  2 
(A7)  a,  =  — —  A,  =  — 

'  '  k  NAt  k  N 


{iVl  fj  cos  <  ^TT-  )  +  4  (-lA, 


N 


N 


j=l 


th 

for  the  k  Fourier  Cosine  Coefficient  of  f.  Thus,  if  we  define 

<A8>  ak  ■  I  f0  +  i  fj  COS  ^ff  +  1  ('1>kfN  • 

j=l 

the  ordinary  Fourier  coefficient  if  given  by 

<A9>  s  ■  £  \ 

and  the  discrete  Fourier  cosine  transform  of  f  is  given  by 
(A10)  Ak  =  At&k  . 
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APPENDIX  B 
PROGRAM  LISTING 


AMPSOE  ••  ANY  PAN'S  POWER  SPECTRAL  DENSITY  ESTIMATOR 


C  R  ARNOLD  12/17/64 


DATA  CARO  1  »  CONTROL 
COL.  1-1(2  *  N 
COL.  11-20  *  K 
COL.  21-30  *  DT 
COL.  31-40  *  SF 


PARAMETERS 

(FIXED  INTEGER  *  NO.  OF  DATA  VALUES  MINUS  ONE) 
(FIXED  INTEGER  =  MAXIMUM  LAG  INDEX) 

(FLT.  PT.  NO.  *  TIME  INCREMENT  BETWEEN  DATA  POINTS) 
( FLT.  PT.  NO.  *  SCALE  FACTOR  FOR  SPECTRAL  DENSITY) 


DATA  CAROS  2-  *  FLT.  PT.  DATA  VALUES,  7  PER  CARD,  COLS.  1-70 


PLOTS  OF  THE  NORMALIZED  AUTOCORRELATION  FUNCTION,  AND  THE  POWER  SPECTRAL 
DENSITY  (ONE  HALF  OF  TWO  SIDED  SPECTRA)  ON  BOTH  A  LINEAR  AND  A  LOG-LOG 
SCALE  ARE  GENERATED  ON  TAPE  A6  WITH  ONE  FILE  FOR  EACH  SET  OF  DATA 


PROGRAM  WILL  TAKE  MULTIPLE  SETS  OF  DATA/EACH  SET  HAVING  A  CONTROL  CARD 


PROGRAM  NORMALLY  TERMINATED  WITH  A  CONTROL  CARD  WITH  A  NEGATIVE 
INTEGER  IN  THE  FIRST  10  COLUMNS 


Cl  MENS  ION  X(  10005 ),R( 1001 ),C( 2001  ),V( 1001 ) , P ( 1001 ), BUFFER (  1024 ) 
CALL  PLOTS(8UFFER( 1024), 1024) 

CALL  PL0T(24. 0,-24. 0,-3) 

CALL  PLOT(0. 0,16. 0,-3) 


READ  IN  CONTROL  PARAMETERS 

1  READ  INPUT  TAPE  2 , 2, NN , KK , DT , SF 

2  FORMAT(2I10,F10.8,F10.5) 

IF ( NN ) 150 , 150,3 

3  WRITE  OUTPUT  TAPE  3,4,NN,KK,DT,SF 

4  FORMAT (1H1, 8X,42HANY  MAN'S  POWER  SPECTRAL  DENSITY  E ST  I MATOR , 3X , 2H* 
1«,3X,10HC  R  ARNOLD, 3X,2H*«,3X, 8 HG ROUP  28, ////9X ,24H INPUT  CONTROL  P 
2ARAMETERS,34X,18HLINC0LN  LABORATORY ,/// /9X , 3HN  =, 15, ///18X,3HK  *tI 
34,///27X,4HCT  *,F8.5, ///36X,4HSF  * , F 10. 5 , ///// ) 


CALCULATE  CONSTANTS 

NP1*NN*1 

NP2*NN*2 

NP3-NNO 

KP1*KK*1 

KP2*KK42 

FK*KK 

FNP1-NP1 

CELTA*0.05 

FMM*0. 43429448 

CP  *2.0«FMM 
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AFPSCE  ••  ANY  FAN  *  S  POWER  SPECTRAL  DENSITY  ESTIMATOR  ••  C  R  ARNOLD 
CF  *3 • 0»FFM 


READ  IN  DATA 

READ  INPUT  TAPE  2 , 6 ,  (  X  (  I  ) ,  I  =  1 »NP1 ) 

6  FORMAT  l 7F  10*7 ) 

WRITE  OUTPUT  TAPE  3,8 

8  FORFAT!1H1,8X,30HFIRST  200  VALUES  OF  INPUT  DATA,//) 
WRITE  OUTPUT  TAPE  3, 10, (XII), 1-1,200) 

10  FORM AT ( 1H  , 8F 15  *6  ) 


CALCULATE  AND  REMOVE  DC  COMPONENT  OF  THE  DATA 

SUF*2 • 0 
CO  12  1*1, NP1 
SUF*SUM4-X  (  I  ) 

12  CONTINUE 
CC*SUM/FNP1 
CO  1 A  1*1, NP1 
xm*x(i)-DC 
14  CONTINUE 

WRITE  CUTPUT  TAPE  3, 16, DC 

16  FORMAT! 1H1,8X,39HTHE  FEAN  OR  DC  COMPONENT  OF  THE  INPUT  *,F9.6,//) 


CALCULATE  AUTO-CORRELATION  FUNCTION 

CO  20  K* 1 , K P 1 
SUF*0.0 
L  I F*NP2-K 
CO  18  1*1, LIM 
KX* I ♦K— 1 

SUF*SUM+X( I )»X(KX) 

18  CONTINUE 
FLIF*LIM 
R!K)*SUM/FLIM 
20  CONTINUE 
V AR*R ( 1 ) 

WRITE  OUTPUT  TAPE  3, 22, VAR 

22  FORMAT! 1H0, //1H0 , // IH0 , 8X , 27HTHE  VARIANCE  OF  THE  INPUT  *,F12.8,//) 

CALCULATE  ANO  PRINT  NORMALIZED  AUTO-CORRELATION 

CO  24  K* 1 , KP 1 
R  <  K ) *R  t  K ) /VAR 
24  CONTINUE 

WRITE  OUTPUT  TAPE  3,26 

26  FORMAT! 1H1, 8X, 31HTHE  NORMALIZED  AUTO-CORRELATION,//) 

WRITE  OUTPUT  TAPE  3, 2 8 , ( R I K ) , K* 1 , KP 1 ) 

28  FORMAT! 1H  .8F15.6) 

GENERATE  TABLE  OF  COSINE  VALUES 


12/17/64 
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AMPSCE  **  ANY  MAN 9 S  POWER  SPECTRAL  DENSITY  ESTIMATOR  ••  C  R  ARNOLD 
KB*2*KP 1 

CNST1=3.14159265/FK 
CO  4£  1*1, KPI 
ALPHA* I— 1 

C(I)*C0SF(ALPHA*CNST1) 

40  CONTINUE 

CO  45  1*1, KK 
KX*KB-I 
C(KX)*C( I ) 

45  CONTINUE 

CALCULATE  COSINE  TRANSFORM 

KX*0 

R  (  1  ) *0. 5*R ( 1  ) 

R(KP1)*0.5*R(KP1) 

MOC*2*KK^l 
50  KB* 1 

CO  55  1*1, KPI 
V(n*C(KB) 

KB*KB4KX 

IF (MOO— KB  >52,55,55 
52  KB*KB*1-M0D 
55  CONTINUE 
SUM*0 • 0 
CO  60  1*1, KPI 
SUM* SUM* V ( I ) *R ( I ) 

6e  CONTINUE 
JX*KX4l 

P( JX)*2,0*DT*SUM 
I  F ( KX— KK 162,65,65 
62  KX*KX*l 
GO  TO  50 

65  R(l)*2.0*R<n 
R(KP1)*2.0*R(KP1 ) 

C  SMOOTH  SPECTRAL  DENSITY  WITH  HANNING  WEIGHTS 

C 

V(  1  )«0.5*P<n*0.5»P  (2) 

CO  66  1*2 , KK 

V(  I)*0.25*P(  I-1K0. 5*P<  I1^0.25«P(  14  1) 

66  CONTINUE 

V ( KP 1 )*0.5*P(KK)40.5«P(KP1) 

C 

C 

C  APPLY  VARIANCE  OR  SCALE  FACTOR  TO  ESTIMATE 


IF  f 

SF ) 67 , 67 , 68 

67 

SF* 

1.0/DT 

GO 

TO  69 

68 

SF* 

SF« VAR 

69 

CO 

70  1*1, KPI 

V  (  I 

)*SF*V( I ) 

P(I 

)*V( I  ) 

70 

CONTINUE 

12/17/64 
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AMPSCE  *«  ANY  PAN'S  POWER  SPECTRAL  OENSITY  ESTIMATOR  *•  C  R  ARNOLD 


PRINT  SPECTRAL  ESTIMATE 
WRITE  OUTPUT  TAPE  3,72 

72  FORMAT! 1H 1 , 8X , 30HTHE  ESTIMATED  SPECTRAL  DENSITY,//) 
WRITE  OUTPUT  TAPE  3 , 28 , ( P ( I) , I  *  1 , KP 1 ) 


GENERATE  TIME  SCALE  FOR  CALCOMP  PLOT  OF  AUTO-CORRELATION 

TMIN*0.0 

SK*KK 

TMAX*SK*DT 
T INC*200. 0*OT 
SK*SK/20.0 
SK  I  P*SK44 • 0 


SCALE  AND  PLOT  NORMALIZED  AUTO-CORRELATION 
R ( KP2 ) *- 1 • 0 

CALL  SCALE(R,KP2, 10. 0,RMIN,OR ) 

CALL  AXIS! 0.0,0. 0, 15HAUTOCORRELATION.15, 1 0. 0 , 90. 0 , RMI N , OR ) 
CALL  AXIS!0.0,5.0,3HLAG,3,SK,0.0,TMIN,TINC) 

13=3 
T*0  •  0 

CO  30  1*1, KPI 
CALL  PLOT ( T , R ( I ) , I  3 ) 

13*2 

T*T*CELTA 
30  CONTINUE 


CALL  PLOT!0. 0,-13. 0,-3) 


GENERATE  FREQUENCY  SCALE  FOR  CALCOMP  PLOT  OF  SPECTRAL  OENSITY 

CM  IN*0.0 
CMAX-0.5/DT 
CQ* 100. 0/TMAX 
C 
C 

C  SCALE  ANO  PLOT  ESTIMATEO  SPECTRAL  OENSITY 

C 

V ( Kp2 ) *0 • 0 

CALL  SCALE! V,KP2, 10.0, PMIN, OP) 

CALL  AXI$<0. 0,0.0, 13HPOWER  OENSITY, 13, 10.0,90.0, PMIN, DP ) 
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AMPSCE  ••  ANY  PAN  *  S  POWER  SPECTRAL  DENSITY  ESTIMATOR  ••  C  R  ARNOLD 

CALL  AXIS(0.0, 0.0, 9H FREQUENCY, 9, SK,0.0»QMIN, DO) 

13*3 
T*2  •  0 

CO  75  I* l *KP l 
CALL  PLOT (T«V(I)vI3) 

13*2 

T*T*CELTA 
75  CONTINUE 
C 
C 

C  CALCULATE  POWER  LEVELS  IN  DB  (IE.  LOG  SCALE) 

C 

co  m  i«i,kpi 

IF ( P ( I  )-0.0e031623)102,103,  103 
102  P(  1 1*0.00031623 
ie3  P(  I)*CP#LOGF(PCI) ) 

V(  I  )*5.0«P( I ) 
ieA  CONTINUE 

WRITE  OUTPUT  TAPE  3,105 
K5  FORMAT!  1H1,8X,18H POWER  LEVELS  IN  DB*///) 

WRITE  OUTPUT  TAPE  3,28, ( VC  I ) , 1*1 ,KP1 ) 

C 

C  CALCULATE  LOG (FREQ )  VALUES 

C 

CF*1.0/(2.0*DT«FK) 

V(1)*0.0 
CO  ie6  1*2, KP1 

v(  n*v(  i-i)4df 
106  CONTINUE 
VC l)*-4.0 
CO  108  1*2, KP1 
V(  I)*CF*LOGF(V(I) > 

108  CONTINUE 
C 
C 

CALL  PLCT(SKIP,0.0,-3) 

C 

C 

CALL  PLOT ( 14.0,0,0,2) 

CALL  PLOT (14.0,12.0,2) 

CALL  PLOT(0. 0,12.0,2) 

CALL  PLOT(0.e,0.0,2) 


CALL  PLOT(4.e,7.0,-3) 

C 

C 

C 

C  PLOT  ESTIMATE  ON  LOG-LOG  PLOT  (BODE  PLOT) 

C 

CALL  SYMBL4(V,P,0.2,lH*»0.0»l ) 

13*3 

CO  110  I*  1 , KP 1 

CALL  PLOT (V(I),P(I),I3) 

110  13*2 
C 
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AMPSCE  «•  ANY  PAN’S  PCJ^ER  SPECTRAL  DENSITY  ESTIMATOR  **  C  R  ARNOLD 


C 

C 


C 

C 

C 


C 

C 

C 

C 


C 

C 

C 

C 


PLOT  AND  LABEL  PCfcER  AXIS 


CALL  AXlS4(0.0,-6.0,  10.0,90.0) 
CALL  SYMBL4 (- 1.20, 3. 9,0.2, 6H420 
CALL  SYMBL4(-1.20, 1.9, 3.2,6H+10 
CALL  SYKBlA(-1.2,-2.1,0.2,6H-10 
CALL  SYMBL4(-l.2,-4.1,0.2,6H-20 
CALL  SYMBL4(-l.2,-6.l,0.2,6H-30 


DB , 0 • 0 , 6 ) 
DB, 0.0,6) 
OB, 0.0, 6) 
DB, 0.0,6) 
06,0.0,6) 


PLOT  ANC  LABEL  FREQUENCY  AXIS 


V(  l)=0.l 
CC  112  1=2,10 
V(  I  )  =V  (  I-D+0.1 
112  CONTINUE 

CC  114  1*11,19 
Ml  I )=V(  I-l  )4l.e 
l  14  CCM  INUE 

CC  116  1*20,20 
V(  I  )  = V (  1-1)410. e 
l  16  CCM  INUE 

CO  118  1=29,37 
V  (  I  )  * V  (  1-1)4100.0 
lie  CCM  INUE 

CC  120  1*1,37 

V(  I  )  =CF  *LCGF ( V (  I ) ) 

120  CCM  INUE 

CALL  PLCT (V(l), -0.1,3) 

CC  130  1*1,36 

CALL  PLCT ( V (  I  ) , 0.0,2 ) 

XXX  =  \M  I  4  l  | 

CALL  PLCT(XXX,0.0,2 ) 

CALL  PLCT ( XXX ,-0 • 1,2) 

130  CCM  INUE 

CALL  SYMBL 4(  40.70,-0. 35, 0.2, AH  1000, 0.0, A) 
CALL  SYMBL4 (45.80, -0.35, 0.2, 3H100, 0.0,3) 
CALL  SYMBL4 ( 4 2 . 90 ,-0 . 3 5 , 0 . 2 , 2H 1 0 , 0 . 0 , 2 ) 
CALL  SYHBLA (-3.23,-0.35,0.2,3H0. 1 ,0.0,3) 


END  PLOT  ANO  GO  BACK  FOR  NEXT  SAMPLE 

1A0  CALL  PLOT(20. 0,6. 0,-3) 

CALL  PLOT(0. 0,0. 0,-3) 

ENC  FILE  6 
GO  TC  1 


END  PROGRAM 

150  ENC  FILE  6 
ENC  FILE  6 
ENC  FILE  6 
CALL  RUN (  1,6) 

CALL  EXIT 

ENC (  l,l,0,0,B,0,l,l,e,0, 0,0,0, 0,0) 
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SAMPLED  LISTED  OUTPUTS 


ANY  MAN’S  POfcER  SPECTRAL  DENSITY  ESTIMATOR  **  C  R  ARNOLD  ** 


INPUT  CONTROL  PARAMETERS 


LINCOLN 


N  =  2e 99 


k  =  2ee 


dt  =  B.eesep 

SF  *  0. 


GROUP  28 


LABORATORY 
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FIRST  200  VALUES  OF  INPUT  OATA 


0.003167 

-0.023012 

0.001821 

-0.015544 

-0.0e5086 

-0.005054 

-0.027427 

-0.017369 

-0.023924 

0.000895 

-0.003641 

0.013763 

0.006324 

0.043818 

0.059757 

0.056577 

0.061506 

0.022936 

0.035334 

0.057114 

0.022063 

0.029074 

0.019606 

-0.023026 

0.002031 


-0.005749 

-0.017448 

0.005840 

-0.012399 

-0.010975 

-0.010342 

-0.026475 

-0.014409 

-0.012422 

-0.002861 

0.000937 

0.011450 

0.006650 

0.040871 

0.053769 

0.053083 

0.070521 

0.021861 

0.038199 

0.046289 

0.023262 

0.033535 

0.011882 

-0.009547 

-0.002342 


-0.011367 

0.001370 

0.001003 

-0.014717 

-0.019002 

-0.010682 

-0.019039 

-0.007360 

-0.006821 

-0.000947 

0.005074 

0.010923 

0.007713 

0.045449 

0.045771 

0.038207 

0.071150 

0.022806 

0.047659 

0.042176 

0.029663 

0.028399 

0.008029 

-0.000501 

-0.009302 
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0.005163 

0.000953 

0.002131 

0.029699 

0.011417 

0.002371 

0.020300 

0.002180 

0.003172 

0.001386 

0.018447 

0.010053 

0.009340 

0.051783 

0.056794 

0.049293 

0.053271 

0.031506 

0.045418 

0.048073 

0*027031 

0.035968 

0.004099 

0.006225 

0.002381 


-0.009799 

-0.002311 

-0.002987 

-0.022186 

-0.012012 

-0.003004 

-0.026474 

-0.004654 

0.004044 

0.005832 

0.012427 

0.008033 

0.012714 

0.039775 

0.047739 

0.052675 

0.050443 

0.033402 

0.043687 

0.041371 

0.032625 

0.022993 

-0.000928 

0.008641 

-0.013952 


-0.003509 

-0.002168 

-0.014390 

-0.027522 

-0.012035 

-0.012269 

-0.024630 

-0.005970 

0.002672 

-0.005780 

0.007912 

0.015801 

0.031227 

0.045987 

0.054025 

0.055173 

0.049653 

0.042664 

0.041664 

0.030891 

0.030824 

0.020388 

0.000658 

0.006158 

-0.011380 


-0.008939 

-0.004476 

-0.014641 

-0.030688 

-0.022996 

-0.013185 

-0.025687 

-0.013058 

0.000501 

-0.006116 

0.004069 

0.011591 

0.036553 

0.041747 

0.059218 

0.056187 

0.042964 

0.041065 

0.042470 

0.029886 

0.035899 

0.021602 

-0.008010 

0.006963 

-0.013439 


-0.019075 

-0.006565 

-0.012506 

-0.022335 

-0.024289 

-0.008765 

-0.034429 

-0.020007 

-0.001307 

-0.002693 

0.013192 

0.015092 

0.042212 

0.050237 

0.050128 

0.058741 

0.038961 

0.038351 

0.047433 

0.028271 

0.038216 

0.012861 

-0.016752 

0.005466 

-0.017056 
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THE  MEAN  OR  DC  COMPONENT  OF  THE  INPUT  =  0.001677 


THE  VARIANCE  OF  THE  INPUT  =  0.002219A6 
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THE  NORMALIZED  AUTO-CORRELATION 


i.eze0e0 

0.988483 

0.977261 

0.905797 

0.892993 

0.880604 

0.811447 

0.801002 

0.790465 

0.734260 

0.724607 

0.713827 

0.644426 

0.633210 

0.621566 

0.555399 

0.545372 

0.535604 

0.481260 

0.472169 

0.463258 

0.413395 

0.405310 

0.397709 

0.351069 

0.344606 

0.338576 

0.305072 

0.299781 

0.293685 

0.260212 

0.255367 

0.250505 

0.217503 

0.211944 

0.206554 

0.173913 

0.168242 

0.162913 

0.140643 

0.136972 

0.133041 

0.111254 

0.108785 

0.107270 

0.100383 

0.097775 

0.095652 

0.082885 

0.082078 

0.080810 

0.073443 

e. 072684 

0.072061 

0.065955 

0.063710 

0.061874 

0.052326 

0.050885 

0.049350 

0.033733 

0.030269 

0.026932 

0.006207 

0.001968 

-0.002705 

-0.024798 

-0.029347 

-0.033794 

-0.058950 

-0.062460 

-0.065813 

-0.089030 

-0.123105 

-0.092410 

-0.096080 

APPENDIX  C 


0.965370 
0.868973 
0.780374 
0.702710 
0.609688 
0.526118 
0.454627 
0.389483 
0.332120 
0.287475 
0.246645 
0.201276 
0.158431 
0. 128752 
0.106156 
0.093491 
0.079162 
0.071322 
0.060037 
0.047386 
0.023246 
-0.007311 
-0.038043 
-0.069140 
-0.099332 


0.954276 

0.857512 

0.770846 

0.691172 

0.597800 

0.516743 

0.445892 

0.381805 

0.326226 

0.281205 

0.241424 

0.196585 

0.154666 

0.125338 

0.104814 

0.091450 

0.078252 

0.070467 

0.058710 

0.044867 

0.019495 

-0.011726 

-0.042253 

-0.072810 

-0.104045 


0.942217 

0.845596 

0.761966 

0.679679 

0.586348 

0.508154 

0.437144 

0.374034 

0.320630 

0.275424 

0.235989 

0.190545 

0.150838 

0.121600 

0.103885 

0.088514 

0.077229 

0.070309 

0.056929 

0.042581 

0.015940 

-0.015393 

-0.046977 

-0.075921 

-0.109162 


0.929966 

0.833682 

0.752618 

0.667647 

0.575479 

0.499441 

0.428739 

0.366265 

0.315612 

0.270487 

0.229746 

0.185269 

0.147311 

0.117824 

0.103440 

0.086273 

0.075797 

0.069132 

0.055543 

0.039943 

0.013107 

-0.018276 

-0.051878 

-0.080768 

-0.113847 


0.918140 
0.822193 
0.7431 73 
0.655792 
0.565433 
0.490440 
0.420930 
0.358213 
0.310338 
0.265387 
0.223451 
0.1 79680 
0.144005 
0.114668 
0.  102465 
0.084227 
0.074900 
0.067951 
0.053550 
0.037485 
0.009692 
-0.021291 
-0.055684 
-0.084992 
-0.118653 


THE  ESTIMATED  SPECTRAL  DENSITY 


93.581613 

72.606593 

36.166934 

1.209590 

1.159778 

1.067970 

0.559497 

0.444835 

0.321608 

0. 105747 

0.075026 

0.125775 

0.073263 

0.121475 

0.136069 

0 .040737 

0.030593 

0.043329 

0.055298 

0.043584 

0.041222 

0.028214 

0.029588 

0.030122 

0.019161 

0.023863 

0.025889 

0.024218 

0.016438 

0.011352 

0.017845 

0.015686 

0.007391 

0.007867 

0.004957 

0.011512 

0.020158 

0.019220 

0.016525 

0.011950 

0.009531 

0.006734 

0.010963 

0.010003 

0.010729 

0.005464 

0.008298 

0.008797 

0.004753 

0.006637 

0.008774 

0.007035 

0.004808 

0.004427 

0.007847 

0.006325 

0.003559 

0.007304 

0.006327 

0.005675 

0.006340 

0.005318 

0.005234 

0.006503 

0.007777 

0.007837 

0.004938 

0.006223 

0.008571 

0.006832 

0.004602 

0.006495 

0.008237 

0.009467 

0.008217 

0.007243 
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15*285943 

0*810181 

0.251954 

0.156942 

0*080033 

0.055089 

0.038913 

0.033064 

0.026414 

0.011369 

0.009969 

0.018343 

0.013963 

0.008498 

0.010580 

0.007196 

0.008139 

0.006850 

0.005693 

0.004520 

0.006487 

0.005883 

0.007431 

0.009681 

0.007272 


7.339891 

0.682259 

0.202903 

0.131297 

0.050445 

0.043175 

0.031075 

0.036628 

0.022617 

0.015252 

0.015455 

0.017637 

0.013575 

0.011610 

0.008943 

0.006039 

0.007234 

0.005963 

0.009216 

0.003878 

0.008283 

0.005000 

0.005986 

0.008133 

0.006735 


4.448726 

0.732379 

0.125660 

0.106723 

0.055309 

0.035928 

0.024761 

0.039529 

0.020893 

0.018363 

0.012838 

0.016475 

0.013100 

0.011404 

0.007991 

0.006110 

0.007532 

0.003903 

0.007620 

0.004805 

0.009058 

0.007896 

0.006882 

0.005146 

0.005779 


2.873741 

0.740237 

0.128743 

0.106586 

0.063480 

0.040190 

0.025874 

0.034146 

0.027652 

0.017022 

0.012544 

0.015812 

0.010667 

0.011458 

0.007483 

0.005953 

0.005802 

0.004759 

0.005949 

0.006551 

0.007462 

0.008132 

0.006119 

0.004287 

0.005465 


1.478441 

0.656483 

0.161164 

0.085253 

0.058536 

0.052992 

0.028505 

0.020532 

0.031167 

0.014108 

0.013433 

0.017177 

0.010726 

0.012473 

0.005267 

0.004892 

0.005917 

0.006591 

0.007101 

0.007397 

0.005628 

0.005715 

0.006530 

0.005610 

0.006471 


PCWER  LEVELS  IN  DB 


19.711905 

18.609760 

15.583116 

0.826383 

0.643749 

0.285589 

-2.522021 

-3.518011 

-4.926736 

-9.757337 

-11.247897 

-9.004056 

-11.351160 

-9.155119 

-8.662407 

-13.900129 

-15.143764 

-13.632233 

-12.572881 

-13.606699 

-13.848740 

-15.495301 

-15.288816 

-15.211140 

-17.175841 

-16.222730 

-15.868859 

-16.158637 

-17.841593 

-19.449128 

-17.484745 

-18.044876 

-21.313200 

-21.041811 

-23.047373 

-19.388446 

-16.955487 

-17.162549 

-17.818516 

-19.226228 

-20.208486 

-21.717168 

-19.600  570 

-19.998612 

-19.694548 

-22.625077 

-20.810076 

-20.556669 

-23.230194 

-21.780471 

-20.568165 

-21.527414 

-23.180190 

-23.539086 

-21.052745 

-21.989533 

-24.486281 

-21.364614 

-21.987853 

-22.459995 

-21.979179 

-22.742128 

-22.811721 

-21.869076 

-21.091760 

-21.058426 

-23.064732 

-22.059978 

-20.669545 

-21.654593 

-23.370324 

-21.873975 
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-21.400802 

-20.238010 

-20.852751 
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11.842922 

8.656896 

6.482357 

4.584476 

1 

-0.914179 

-1.660510 

-1.352640 

-1.306292 

-1 

-5.986784 

-6.927122 

-9.008046 

-8.902752 

-7 

-8.042615 

-8.817438 

-9.717423 

-9.723003 

-10 

10.967301 

-12.971856 

-12.572010 

-11.973599 

-12 

12.589360 

-13.647670 

-14.445651 

-13.958782 

-12 

14.099100 

-15.075921 

-16.062290 

-15.871309 

-16 

14.806384 

-14.361810 

-14.030878 

-14.666584 

-16 

15.781597 

-16.455670 

-16.800065 

-15.582720 

-15 

19.442949 

-18.166800 

-17.360509 

-17.689894 

-18 

20.013465 

-18.109404 

-18.915038 

-19.015637 

-18 

17.365269 

-17.535735 

-17.831777 

-18.010159 

-17 

18.550284 

-18.672472 

-18.827414 

-19.719727 

-19 

20.706888 

-19.351708 

-19.429600 

-19.408992 

-19 

19.755076 

-20.485268 

-20.974099 

-21.259058 

-22 

21.429288 

-22.190133 

-22.139465 

-22.252278 

-23 

20.894531 

-21.406152 

-21.230669 

-22.364450 

-22 

21.642945 

-22.245610 

-24.085538 

-23.225269 

-21 

22.446638 

-20.354578 

-21.180660 

-22.255781 

-21 

23.448793 

-24.113971 

-23.182712 

-21.837173 

-21 

21.879473 

-20.818088 

-20.429798 
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-22 

22.303816 

-23.010666 

-21.026017 

-20.897881 

-22 

21.289626 

-22.228683 

-21.622881 

-22.133362 

-21 

20.140699 

-20.897345 

-22.885211 

-23.678906 

-22 

21.383698 

-21.716848 

-22.381790 

-22.623821 

-21 

698040 

827768 

927331 

692918 

325783 

757925 

450785 

875745 

063045 

505206 

718239 

650464 

695793 

040171 

784139 

105397 

279319 

810688 

486724 

309591 

496240 

429556 

850578 

510258 

890021 
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